How to use wfc3_dash on DASH data

Learning Goals:

By the end of this tutorial, you will:

  • Create difference files, association tables, and segmentation maps using wfc3_dash.
  • Subtract background and fix cosmic rays from newly generated FLTs.
  • Align reads to each other for a final product.

Table of Contents:

Introduction
1. Imports
2. Download relevant data
3. Run DASH

4. Plot original IMA vs. DASH pipeline science result
5. Conclusions
Additional Resources
About the Notebook
Citations

Introduction

This notebook is the first in a new Drift And SHift (DASH) pipeline workflow developed to ease the process of reducing DASH data. The pipeline is customizable, able to be changed according to scientific goals of the user, and this first tutorial walks the user from data download to a finished product ready for science analysis. For more information, see Momcheva et. al 2016 and WFC3 ISR 2021-01.

1. Imports

This notebook assumes you have created the virtual environment in WFC3 Library's installation instructions.

We import:

  • os for setting environment variables
  • glob for querying through directories
  • numpy for handling array functions
  • matplotlib.pyplot for plotting data
  • astropy for astronomy related functions
  • astroquery.Observations for downloading data
  • drizzlepac.astrodrizzle for combining images
  • reduce_dash for reducing DASH data
In [ ]:
%matplotlib notebook

import os
from glob import glob

import numpy as np
import matplotlib.pyplot as plt 

from astropy.io import fits 
from astropy.table import Table
from astropy.io import ascii
from astroquery.mast import Observations

from drizzlepac import astrodrizzle 

from reduce_dash import DashData

2. Download relevant data

Retrieve the table of observations associated with 15238.

In [ ]:
obsTable = Observations.query_criteria(proposal_id=['15238'], obs_id=['IDNM0J030'])

Get the full list of products associated to the table and restrict the list to IMA files.

In [3]:
product_list = Observations.get_product_list(obsTable)
BM = (product_list['productSubGroupDescription']  == 'IMA') 
product_list = product_list[BM]

product_list.show_in_notebook(display_length=5)
Out[3]:
Table length=2
idxobsIDobs_collectiondataproduct_typeobs_iddescriptiontypedataURIproductTypeproductGroupDescriptionproductSubGroupDescriptionproductDocumentationURLprojectprvversionproposal_idproductFilenamesizeparent_obsiddataRightscalib_level
026144024HSTimageidnm0jodqDADS IMA file - Intermediate Mult-Accum WFC3/NICMOSSmast:HST/product/idnm0jodq_ima.fitsAUXILIARY--IMA--CALWF33.5.2(Aug-12-2020)15238idnm0jodq_ima.fits16826112026144060PUBLIC2
126144027HSTimageidnm0jofqDADS IMA file - Intermediate Mult-Accum WFC3/NICMOSSmast:HST/product/idnm0jofq_ima.fitsAUXILIARY--IMA--CALWF33.5.2(Aug-12-2020)15238idnm0jofq_ima.fits16826112026144060PUBLIC2

Choose a single exposure file to work on. In this example, we choose the first exposure. To create usable data, you will have to follow this work flow on all individual IMA files in your dataset.

In [ ]:
myID = product_list['obsID'][0:1]

Download the IMA and FLT files for that exposure. The standard pipeline-FLT will be used for comparison with the detrended final product.

In [ ]:
download = Observations.download_products(myID,mrp_only=False,productSubGroupDescription=['IMA','FLT'])

Display the results of the download operation.

In [6]:
download
Out[6]:
Table length=3
Local PathStatusMessageURL
str99str8objectobject
./mastDownload/HST/hst_15238_0j_wfc3_ir_f139m_idnm0jod/hst_15238_0j_wfc3_ir_f139m_idnm0jod_flt.fitsCOMPLETENoneNone
./mastDownload/HST/idnm0jodq/idnm0jodq_ima.fitsCOMPLETENoneNone
./mastDownload/HST/idnm0jodq/idnm0jodq_flt.fitsCOMPLETENoneNone

Read the files that were just downloaded locally. In addition, have the path be just the rootname, i.e. without the file extension.

In [ ]:
localpathtofile = download['Local Path'][1][:-8]
localpathtofile

original_ima = fits.open(localpathtofile+'ima.fits')
original_flt = fits.open(localpathtofile+'flt.fits')
original_ima.info()

Print the number of samples and plot the individual reads of the IMA file.

Note: the individual 'SCI' extensions are stored in reverse order, with 'SCI', 1 corresponding to the last read.

In [8]:
nsamp = original_ima[0].header['NSAMP']
print('NSAMP',nsamp)

fig,axarr = plt.subplots((nsamp+3)//4,4, figsize=(9,3*((nsamp+3)//4)))

for i in range(1,4*((nsamp+3)//4)+1):

    row = (i-1)//4
    col = (i-1)%4
    if (i <= nsamp):
        immed = np.nanmedian(original_ima['SCI',i].data)
        stdev = np.nanstd(original_ima['SCI',i].data)
        axarr[row,col].imshow(original_ima['SCI',i].data,clim=[immed-.3*stdev,immed+.5*stdev],cmap='Greys',origin='lower')
        axarr[row,col].set_title('SCI '+str(i))
        axarr[row,col].set_xticks([]) 
        axarr[row,col].set_yticks([]) 
    else:
        fig.delaxes(axarr[row,col])

fig.tight_layout()

NSAMP 16

3. Query CRDS for reference files

Before running reduce_dash, we need to set some environment variables for several subsequent calibration tasks.

We will point to a subdirectory called crds_cache/ using the IREF environment variable. The IREF variable is used for WFC3 reference files. Other instruments use other variables, e.g., JREF for ACS.

In [ ]:
os.environ['CRDS_SERVER_URL'] = 'https://hst-crds.stsci.edu'
os.environ['CRDS_SERVER'] = 'https://hst-crds.stsci.edu'
os.environ['CRDS_PATH'] = './crds_cache'
os.environ['iref'] = 'iref/'
if not os.path.exists('iref'):
    os.mkdir('iref')

The code block below will query CRDS for the best reference files currently available for these datasets and update the header keywords to point to these new files. We will use the Python package os to run terminal commands. In the terminal, the line would be:

crds bestrefs --files [filename] --sync-references=1 --update-bestrefs

...where 'filename' is the name of your fits file.

In [ ]:
ima_files = glob('*_ima.fits') 

for file in ima_files:
    command_line_input = 'crds bestrefs --files {:} --sync-references=1 --update-bestrefs'.format(file)
    os.system(command_line_input)

3. Run DASH

Run the DASH pipeline for a single exposure. This procedure showcases the capabilities and customization options of the DASH pipeline.

Note: the following will only work if you are using the notebooks inside of the Notebook directory. wfc3_dash submodule will be properly packaged and installed within the wfc3_tools module sometime in the future.

If you move the notebooks and want to use them elsewhere, you can still provide a temp_path to the dash codes and remove the comments below.

In [ ]:
#import sys

#module_path = os.path.abspath(os.path.join('..'))
#if module_path not in sys.path:
#    sys.path.append(module_path)

#tmp_path = ".../wfc3_dash/wfc3_dash"
#if tmp_path not in sys.path:
#    sys.path.append(tmp_path)

3.1 Create DashData object

We use the both IMA and FLT extensions of our local image to create a DashData object.

In [ ]:
myDash = DashData(localpathtofile+'ima.fits', flt_file_name=localpathtofile+'flt.fits')
print(myDash.root)

3.2 Create difference files from reads

A difference (diff) file contains the counts accumulated between two reads of the IMA file. The diff files are written to disk in a directory named ./diff under the current working directory (cwd). In creating diff files from the readouts of the IMA, the first difference, between the 1-st and 0-th read is ignored because of its very short exposure time of 2.9 seconds, resulting in a noisy image.

In order to create the best possible results, the split_ima() method uses the bestrefs function from CRDS to ensure all reference files are up to date and available.

In [ ]:
myDash.split_ima()

Print the number of diff files and plot the diff files.

In [14]:
ndiff = len(myDash.diff_files_list)
print('Number of diff files',ndiff)

if ndiff > 4: 
    fig,axarr = plt.subplots((ndiff+3)//4,4, figsize=(9,3*((ndiff+3)//4)))

    for i in range(4*((ndiff+3)//4)):

        row = (i)//4
        col = (i)%4
        if (i < ndiff):
            diff_i = fits.open(myDash.diff_files_list[i]+'_diff.fits')
            immed = np.nanmedian(diff_i['SCI'].data)
            stdev = np.nanstd(diff_i['SCI'].data)
            axarr[row,col].imshow(diff_i['SCI'].data,clim=[immed-.3*stdev,immed+.5*stdev],cmap='Greys',origin='lower')
            axarr[row,col].set_title('Diff:'+str(i+1))
            axarr[row,col].set_xticks([]) 
            axarr[row,col].set_yticks([]) 
        else:
            fig.delaxes(axarr[row,col])
else:
    fig,axarr = plt.subplots(1,ndiff,figsize=(15,15))
    for i in range(ndiff):
        immed = np.nanmedian(diff_i['SCI'].data)
        stdev = np.nanstd(diff_i['SCI'].data)
        diff_i = fits.open(myDash.diff_files_list[i]+'_diff.fits')
        axarr[i].imshow(diff_i['SCI'].data,clim=[immed-.3*stdev,immed+.5*stdev],cmap='Greys',origin='lower')
        axarr[i].set_title('Diff:'+str(i+1))
        axarr[i].set_xticks([]) 
        axarr[i].set_yticks([]) 

fig.tight_layout()

Number of diff files 14

3.3 Create an association table

This file mimics a typical association file for dithered exposures, which is used by AstroDrizzle to align and stack multiple exposures taken at the same sky position with small dithers.

We exploit the fact that a WFC3/IR exposure taken under gyro control can be effectively split into individual pseudo-exposures (the diff images we created in Section 3.2). From there, AstroDrizzle can treat such pseudo-expsoures as individual dithers, and combine them into a single exposure.

In [ ]:
myDash.make_pointing_asn()

Show the content of the asn file, which was created in ./diff.

In [ ]:
asn_filename = 'diff/{}_asn.fits'.format(myDash.root)
asn_table = Table(fits.getdata(asn_filename, ext=1))
asn_table.show_in_notebook()

3.4 Create a segmentation map

Create segmentation map from original FLT image to assist with background subtraction and fixing of cosmic ray flags using create_seg_map(). This method makes a directory called ./segmentation_maps, which holds the outputs.

In [ ]:
myDash.create_seg_map()

Plot segmentation map.

In [18]:
rootname = myDash.root
segmap_name = ('segmentation_maps/'+ rootname + '_seg.fits')
segmap = fits.getdata(segmap_name)

fig = plt.figure(figsize=(6, 8))
plt.title(segmap_name)
plt.imshow(segmap, origin='lower', vmin=0, vmax=1, cmap='Greys_r')
Out[18]:
<matplotlib.image.AxesImage at 0x7ffd550011d0>

Print and read source list.

In [ ]:
sourcelist_name = ('segmentation_maps/' + rootname + '_source_list.dat')
sourcelist = ascii.read(sourcelist_name)
print(sourcelist)

Let's create a segmentation map and source list from the difference files. We need to make source lists from our difference files created from the IMA so that TweakReg can better align these difference files to catalogs, each other, etc.

First, generate a list of difference files that contain the full path name.

In [ ]:
diffpath = os.path.dirname(os.path.abspath('diff/{}_*_diff.fits'.format(rootname)))
cat_images = sorted([os.path.basename(x) for x in glob('diff/{}_*_diff.fits'.format(rootname))])
sc_diff_files = [diffpath + '/' + s for s in cat_images]

Then, create a difference segmentation map using diff_seg_map() and the diff files.

In [ ]:
myDash.diff_seg_map(cat_images=sc_diff_files)

Plot the segmentation map from a diffrence files.

In [22]:
segmap_name = ('segmentation_maps/' + rootname + '_01_diff_seg.fits')
segmap = fits.getdata(segmap_name)
fig = plt.figure(figsize=(6, 8))
plt.title(segmap_name)
plt.imshow(segmap, origin='lower', vmin=0, vmax=1, cmap='Greys_r')
Out[22]:
<matplotlib.image.AxesImage at 0x7ffd4243ea10>

3.5 Subtract background from difference files

Subtract background from the individual reads taken from the original IMA file using the DRZ and SEG images produced in the background subtraction of the original FLT.

By default, subtract_background_reads() will subtract the background and write it to the header. By setting subtract=False, the background will not be subtracted and will only be written to the header. In addition, setting reset_stars_dq=True will reset cosmic rays within objects to 0 since the centers of the stars are flagged.

In [ ]:
myDash.subtract_background_reads()

3.6 Fix cosmic rays

Now, we can use fix_cosmic_rays() to reset cosmic rays within the segmentation maps of objects and use L.A.Cosmic to find them again.

In [ ]:
myDash.fix_cosmic_rays()

3.7 Align reads to each other

Align reads from the IMA to one another by aligning each difference file to the first diff file.

Listed below are all the parameters available to myDash.align(). align() uses TweakReg to update the WCS information in the headers of the diff files, then drizzles the images together using AstroDrizzle. There are more parameters available to users when working with TweakReg and AstroDrizzle that could be an integral part of the workflow for users of DASH. The example below lists the default values set for every input:

myDash.align(self, subtract_background = True, align_method = None, ref_catalog = None, create_diff_source_lists = True, updatehdr = True, updatewcs = True, wcsname = 'DASH', threshold = 50., cw = 3.5, searchrad = 20., astrodriz = True, cat_file = 'catalogs/diff_catfile.cat', drz_output = None, move_files = False)

Refer to documentation to customize parameters for TweakReg and AstroDrizzle.

Note: the error UnboundLocalError: local variable 'sig' referenced before assignment can be solved by lowering threshold parameter.

In [ ]:
myDash.align(updatehdr=False, updatewcs=True, astrodriz=False)

Print the shifts file to analyze how well the alignment went. Do not update header until shifts, as seen in the xrms and yrms columns, are satisfactory. Further information about the outputs in the shift file and what is 'satisfactory' can be found in the Drizzlepac Handbook.

In [ ]:
shift_file = glob('shifts/shifts_*.txt')
shift_file_name = shift_file[0]


shift_tab = Table.read(shift_file_name,
                       format='ascii.no_header',
                       names=['file','dx','dy','rot','scale','xrms','yrms'])

formats = ['.2f', '.2f', '.3f', '.5f', '.2f', '.2f']
for i, col in enumerate(shift_tab.colnames[1:]):
    shift_tab[col].format = formats[i]
shift_tab

Let's align our images with a threshold of 20, and update the headers and WCS information.

In [ ]:
myDash.align(threshold = 20.)

4. Plot original IMA and DASH pipeline science result

Plot the final DRZ image and compare to the original IMA.

In [28]:
sci_name = myDash.root + '_drz_sci.fits'
og_flt_name = 'mastDownload/HST/' + myDash.root + '/' + myDash.root + '_ima.fits'
sci = fits.getdata(sci_name)
og_flt = fits.getdata(og_flt_name)

fig = plt.figure(figsize=(9, 4))
ax1 = fig.add_subplot(1,2,2)
ax2 = fig.add_subplot(1,2,1)

ax1.set_title('DASH Pipeline Reduced Science File')
ax2.set_title('Original IMA (not reduced using pipeline)')

ax1.set_xlim(-10,1120)
ax2.set_xlim(-10,1120)

ax1.set_ylim(-10,1050)
ax2.set_ylim(-10,1050)

ax1.imshow(sci, vmin=0, vmax=40, cmap='Greys_r', origin='lower')
ax2.imshow(og_flt, vmin=0, vmax=40, cmap='Greys_r', origin='lower')
plt.show()

5. Conclusions

Thank you for walking through this notebook. Now using WFC3 data, you should be more familiar with:

  • Creating difference files, association tables, and segmentation maps using wfc3_dash.
  • Subtracting background and fixing cosmic rays from newly generated FLTs.
  • Aligning reads to each other for a final product.

Congratulations, you have completed the notebook!

Additional Resources

Below are some additional resources that may be helpful. Please send any questions through the HST Helpdesk.

About this Notebook

Authors: Catherine Martlin; WFC3 Instrument Team

Updated on: 2021-10-07

Citations

If you use numpy, matplotlib, astropy, astroquery, or drizzlepac for published research, please cite the authors. Follow these links for more information about citing the libraries below:

Top of Page Space Telescope Logo

In [ ]: